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Ptychography is a popular technique to achieve diffraction limited resolution images of a two or 
three dimensional sample using high frame rate detectors. We introduce a relaxation of common 
projection algorithms to account for instabilities given by intensity and background fluctuations, 
position errors, or poor calibration using multiplexing illumination. This relaxation introduces 
an additional phasing optimization at every step that enhances the convergence rate of common 
projection algorithms. Numerical tests exhibit the exact recovery of the object and the noise when 
there is high redundancy in the data. 



I. INTRODUCTION 

In a ptychography experiment 1 6 (see Fig. 1), a two 
dimensional small beam with distribution w(r) of dimen- 
sion m x m illuminates an unknown object of interest 
ijj{r + x^)), where x^) is used to denote the location of 
the subregion. For simplicity we consider square matri- 
ces, generalization to non-square matrices is straightfor- 
ward but requires more indices and complicates notation. 
We assume that a sequence of k diffraction images a 2 ^ (q) 
are collected as the position x^) of the object is rastered. 
Each frame contains m x m pixels. We consider the 
case in which is the magnitude of the discrete two 
dimensional discrete Fourier transform (DFT) of the un- 
known object. If we let J 7 be a 2D DFT operator, the 
relationship among a^), w(r) and -0(r + x^)) can be ex- 
pressed as follows: 



<Hi)(q) = ^(r)^(r + x (i) ) 
r = rm, q = 

-F/ = y> iqr /(r), 



(1) 



m = (/i, v) , /i, v = (0 . . . m — 1), 

where r is a lengthscale, and the sum over r is given on 
all the indices mxmofr. 

As x is rastered on a typically coarser grid, r -fx spans 
a finer grid of dimension n x n, where n > m. We denote 
with subindex i = (1 • • • k) the sequence of raster points 

X(<). 

We introduce the illumination operator associated 

with translation that extracts a frame out of ^, and 
scales the frame point-wise by the illumination function 
w(r): (see Fig. 1) 

%)(r) = iu(r)^(r + x (i) ) = Q (i )(r)^, 



where represents the frames extracted from ijj and 
multiplied by the probe w(r). 

Using linear algebra notation, we represent ^ as a vec- 
tor of length n 2 , tj; G C n The moving beam associ- 
ated with the illumination function w(r) can be repre- 
sented as an m 2 x n 2 sparse "illumination matrix" 
(i = l,2,-.. ,fe). 

The relationship between the diffraction measurements 
collected in a ptychography experiment and the unknown 
object to be recovered can be represented compactly as: 



1^1, 



(2) 



if we stack the diffraction measurements into a long 
vector a, and define various matrices as follows: 
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The objective of the ptychographic reconstruction 
problem is to find tj; given a from (Eq. 2). This is of- 
ten formulated using a "divide and conquer" approach 
referred to as projection algorithms, iterative transform 
methods, or alternating direction methods 7 . One formu- 
lates the relationship (Eq. 2) as: 

a=\Fz\, (3) 

z = QV>. (4) 

These algorithms are often defined in terms of two projec- 
tion operators Pp and Pq that project onto the solution 
z to Eqs. 3 and 4 that is closest to the current estimate 
described in Section II. 

Alternative approaches include formulating the prob- 
lem as: 

min||a- |FQ^|| 
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FIG. 1: (Left) Experimental geometry in ptychography: An unknown sample with transmission ip(r) is rastered 
through an illuminating beam w(r), and a sequence of diffraction measurements \a^\ 2 is recorded on an area 
detector as the sample is rastered around. The point- wise product between illuminating function and sample 

z (i)( r ) — ^(r)^ 1 * + x («)) * s re l a t e d to the measurement by a Fourier magnitude relationship a^)(q) = LFz^) . 



and solve by standard unconstrained minimization al- 
gorithms such as projected gradient, conjugate gradi- 
ent, Newton and quasi- Newton methods 8 ' 9 . Another 
approach proposes a convex relaxation of the quadratic 
problem by lifting to an n 2 x n 2 space 10 12 . 

In the following section II we will describe the stan- 
dard operators commonly used in the literature. In sec- 
tion III we will introduce an intermediate variable c^\, 
to account for intensity fluctuations, replacing Eq. (4) 
with = Q(i)ip, z = (1, . . . , A;). In section IV we ex- 

tend this approach to situations in which intensity data 
is added incoherently or multiplexed. In section V we 
extend this approach to situations in which the data is 
corrupted by position errors. In Section VI we introduce 
an unknown background b(q), different for every pixel, 
but constant throughout the acquisition (or different 
for every frame z, but constant throughout each frame). 
In Section VII we report the following numerical results: 

• Exact reconstruction with intensity fluctuation 
given by the coefficients (see Fig. 4). 

• Accelerated convergence (Fig. 5) even when no in- 
tensity fluctuation is present in the data by using 
only the phase of (Table II) 

• Exact reconstruction with multiplexing using 4 si- 



multaneous illuminations adding incoherently on 
the detector, with perturbation of the amplitudes 
(Fig. 6) 

• Exact reconstruction of fluctuating background 
noise independent from the sample (Figs. 8,9). 

• Position recovery (Figs. 7) of the illuminating 
probe using iterative Taylor expansion optimiza- 
tion. 



II. STANDARD PROJECTION ALGORITHMS 

The projection operator Pp mentioned in the previous 
section is often known as the Fourier magnitude projec- 
tion operator. Applying this operator to a vector z yields 



Fz 

\Fz\ 



• a. 



(5) 



where division and product are intended as element-wise 
operations. It is easy to verify that \Ppz\ = \a\ and there- 
fore Ppz satisfies Eq. (3) for any z. Pp is a projection 
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in the sense that 

PpZ = argm_in||z (i) -z (i) ||, 
subject to \Fz\ = a, 



(6) 



where || || denotes the Euclidean norm. The matrix Q de- 
fines an orthogonal projection operator Pq that projects 
any vector in C km2 onto the range of Q: 



P Q = Q(Q*Q)- 1 Q* ) 



(7) 



In the simple alternating projection algorithm, the ap- 
proximation to the solutions of (4) and (7) are updated 
by: 

z (/+1) = [PqPf]z (/) 

^+1) = (Q*Q)-lQ* z (*+l) B 

Here ip^ is the running estimate of ?p that solves the 
linear subproblem: 

^ = argmin||*^ -Q^||, (8) 

One may also introduce a regularization factor e and up- 
date the running estimate as: 

^ = (Q*Q + e^Q* (z( £+1 ) + c^O) . 

with typically = 0. The projection operator Pq can 
be expressed in terms of the solution to Eq. 8: 

Pqz^ = Q(Q*Q)- 1 Q*z^ = Q^ } . 

(Q*Q) _1 is a normalization factor corresponding to the 
sum of all the intensities of the illumination function. Q* 
is the operator that scales all the frames by the con- 
jugate of the probe w*(r) and translates them by — x^) 
onto the image ip: 



A number of heuristic algorithms has been 
proposed 13 ' 14 , a few examples are given in Tab. I, 
with j3 G [0, 1] is a relaxation parameter. Very recently, 
an alternating direction method (ADM) designed to 
work with a special augmented Lagrangian function 7 . 
This function is minimized by applying a block coordi- 
nate descent scheme (or alternating search directions) 
akin to these projection operators. 



III. FLUCTUATING INTENSITIES, AND 
AUGMENTED PHASE RETRIEVAL 

One of the practical issues one may face in ptychogra- 
phy is the intensity fluctuation among different diffrac- 
tion frames introduced by instabilities in the light source, 
optics and shutters. Such fluctuation can be accounted 



for by introducing a scaling factor q G C for each diffrac- 
tion frame. By allowing Ci to be complex, we introduce 
an optimization of the relative phase between frames. As 
a result, the definition Z{ is modified so that the equation 



(9) 



holds for i = 1,2, The phase of the scaling factor 

Ci is a relaxation parameter that enables us to improve 
the convergence of iterative algorithms. As we will see 
below, the phase of c enables us to find the intersection 
between Pp and Pq more easily since a constant phase 
factor multiplying each frame is insensitive to the Fourier 
magnitude projection. By optimizing the phase among 
each frame, we solve the phase at a resolution given by 
the step size between frames. 

Since both q and tj; are unknown in (9), the solution to 
(9) is clearly not unique. To exclude the trivial solution 
Ci = 0, for i = 1,2,...& and ip = 0, we introduce an 
additional constraint and solve 



(VWi, c (i) ) = arg min V \Q^ip - c 



(10) 



subject to 



We can also eliminate ip and solve the following problem: 
arg min - P Q )Diag (c) z|| 2 + A a - fc) 

Since it is sufficient to determine the relative difference 
among different q's, the choice of 7 is arbitrary. A rea- 
sonable choice is 7 = k. We can show (see the Appendix) 
that 



(ii) 



where a is chosen to ensure that J2i=i c i = k holds, and 
the matrix H is given by: 



( C1 \ 




/1\ 


c = aH\l, c= : j 


1 = 








w 



H i,j 
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( P Q)i,j = Q'Q Q - 



1 :Q* 



(12) 



Once c is available, the least squares approximation to 
ip can be written as 



^ = (Q*Q) _1 Q*^, D c 



(c x I 



To modify projection algorithms to account for inten- 
sity fluctuations, we replace the operator Pq used in the 
standard projection algorithms listed in Table I by an 
augmented projection operator Pq defined as: 

Po = D^PqD*. 
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projection algorithm 


updating formula z^ +1 ^ = 


Alternating Projection 15 


[PqPf] z (<) 


HIO 15 


[PqPf + (I-Pq)(I-PPf)]x w 


Difference Map 5 


[PfPq + (I - Pf)(I - PPq)]z W 


RAAR 14 


\2pP Q P F + (1 - 2p)P F + p{P Q - I)] z<" 



TABLE I: Popular fix-point algorithms used in phase retrieval 



An alternative modification is to redefine the probe as 
Q c = D~ 1 Q by incorporating the scaling factors q, and 
write: 

Pq = q c (q:q c ) _1 q: 

Although the construction of Pq is motivated by the 
need to account for intensity fluctuations among different 
diffraction frames in the measured data, it turns out to 
be a useful technique for accelerating the convergence of 
projection algorithms even when no intensity fluctuation 
is present in the data. This is because we perform an 
additional phasing optimization around the Pq projec- 
tion that is invisible to the Fourier magnitude projection. 
This enables the intersection between the constraints as- 
sociated with these Pp and Pq to be detected more eas- 

iiy. 

IV. MULTIPLEXING AND INCOHERENT 
MEASUREMENTS 

A practical issue in ptychography are the strict require- 
ments of the experimental geometry to achieve high qual- 
ity data. For example the need for stable, coherent illu- 
mination of the sample, careful control of vibrations, and 
limited detector response function all contribute to limit 
the specifications of a ptychographic microscope. Here 
we consider the case of multiplexing incoherent measure- 
ments to relax some of the experimental constraint. The 
sum of intensities, rather than amplitudes is a result of 
the incoherent addition with varying experimental pa- 
rameters during the measurement. By relaxing the model 
of the measurement process and solving the optimization 
problem described in the previous section, we attempt 
to recover unknown instantaneous fluctuations of the pa- 
rameters and improve convergence. 

The measurement model is as follows: 



where Q, is a sparse matrix that represents an averag- 
ing operation of the intensities. For example, a single 
exposure may represents the sum of the intensities gen- 
erated by a an illumination beam that translates during 
the exposure, or may represent a binned version of the 
high dimensional signal. We consider Z{ the highly re- 
dundant set of frames generated for all the positions of 



the illumination function during an exposure. More ef- 
ficient parameterization of the data are also applicable, 
(see next section). The operator Q then simply adds the 
intensity of the frames associated with an exposure. 

The projection operator associated with this problem 
can be written as follows: 

/ a 2 

where ^* is an operator that copies the entries over all 
the frames that contribute to an exposure. Replacing 
this operator in the reconstruction process is subject of 
recent interest by several groups 16 . 

However, since fluctuations about an average are un- 
constrained by the measurement, Pp is insensitive to 
modifications of the relative amplitude of the frames. We 
can achieve higher performance if we optimize the the rel- 
ative amplitude and phases of the frames as described in 
the previous section. 

Numerical tests described in section VII show the ex- 
act recovery (within numerical precision) of the object 
and a multiplexing array of beam positions averaged in- 
coherently with errors in the calibration of the amplitude 
factors. The ability to recover the relative amplitude of 
a redundant set of frames that are averaged during the 
measurement enables us also to identify which instance of 
the experimental parameters occurred. The augmented 
projection described in the previous section attempts to 
do just want we need to recover or calibrate the ampli- 
tude coefficients of our multiplexing array. 

V. POSITION ERRORS 

We consider the case in which the probe w is trans- 
lated from the input coordinate by an unknown distance 
£. We call the unknown illumination matrix used to 
generate the data. To determine the illumination ma- 
trix, we determine the parameter £ so that the error sq^ 
is minimized: 

argmin || [I - PqJ z|| 2 (14) 

We restrict initially to 1-dimensional translations. Given 
the illumination function w, we can compute the first 
and second order derivatives with respect to transla- 
tion. We denote by Qi,Ri,Si the illumination matrices 
that extract a frame out of an image and multiplies by 
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(u>(x), d x w(x), d 2 w{-x)) respectively. Q, R, S are tall and 
skinny matrices of the same size as Q discussed earlier, 
with identical location of the non-zero entries. The probe 
is perturbed to second order as follows: 

Q^ = Q + £R + £ 2 S. 

where £ is diagonal and real, and represents the trans- 
lation distance different for every frame i. When the 
translation £ has two components, we write: 



M:;) R =(*:) s =( 



where R2, • • • £22 are the illumination matrices that 
extract a frame out of an image and multiply by (w(x), 
d Xl w(x), d X2 w(x.), . . . , d 2 2 w(-x)) respectively. The second 
order term in the Taylor expansion becomes the sum over 
all dimensions: 



C 2 s 

Sx 



\ ( S 12 



Setting d|J • || and d\ 



€2 1 



-S21). 

I to gives (see appendix B): 




, (15) 



using the definition zri j ... j s22 



= [Rl, . . . , S22] qTqQ*Z, 

z = [/ — Pq]z, and where the matrices #1,2, x are defined 
as: 

( ff i)»7 = ( z Rh z Rh - 2zJz Sll4 ) S i:i - z* (O n ) Zj + cc, 
( Z R 2 , Z R 2 * - 2z -* z s 22i ) *ij - z* (022)^ Zj + cc, 
( z Rh z R 2 , - 2zJz 5x J ^ - z* (Ox)^- zj + cc, 



where 



(0 22 )„ = (R 2 ) i «4q (R2) 



Q*Q v— i/j 

)t Q*Q V- -v:/ , 

(O x ) i ^(R 1 )^(R 2 )*. 

The system of equations can be solved efficiently 
by sparse linear algebra solvers. The entries of the 
equation are given by the dot product between frames 
(z, z, zr x , . . . , zs 22 ) with partial overlap and scaling fac- 
tors given by RJ j2 qtqR-1,2- The terms zzs are higher 
order corrections close to the solution and can be ne- 
glected in practice. In Section VII we will show that this 
method can recover the position perturbations to numer- 
ical accuracy when the perturbations are a fraction than 
the probe width. 

VI. BACKGROUND NOISE 

Here we consider an unknown offset b(q) (background) 
different for every pixel but constant throughout the ac- 
quisition: 



(16) 



A trivial variation of the problem is when b^(q) is dif- 
ferent for every frame but the same for every pixel q). 

At each iteration, we solve the following offset mini- 
mization problem with an additional scaling parameter: 



min 



Z(i)\ 



z {i) 



subject to \z {i) \ 2 = r] (a 2 {i) - b^j , 



(17) 



where 77 G M m is a shrinkage parameter that accounts for 
the fact that |i^| 2 is on average smaller than a 2 . This 
is because z^ is obtained from a sequence of linear pro- 
jections that reduce the overall norm. Since z is smaller, 
the solution to the off-set projection problem 17 is bi- 
ased towards a larger offset. Introducing the shrinkage 
parameter 77 equal for every frame provides the flexibility 
to avoid this problem. 

By solving for 77 first, we obtain the first and second 
order terms: 



E< d 



z {i) 



(18) 



(i) 



where = a 2 ^ — b. Solving for b for a fixed 77 gives 



6 «)_ & (*-D = l^ d(i) _ 



z (i) 



2 1 



(4) 



z {e) 
z (i) 

z (i) 



di 



z (t) 
z (i) 



To avoid strong perturbations, however, we set ?y(q) = .8 
if ?y(q) < 0.8. When optimizing for a fluctuating offset 
(b(i)(<l) = fr(i)(0) constant for every frame), we simply 
replace the sum over i with the sum over q. The update 
of z is then computed as a regular Fourier magnitude 
projection operator with an intensity offset: 



p(«(«)-» (0 ) ; W_ ; W 



I ( l ) I 



where we used the notation P = TPT* . 

In the following section we will show that common pro- 
jection methods can recover the background even if the 
SNR is much smaller than 1. 



VII. NUMERICAL TESTS 

The object used to simulate the diffraction pattern is 
obtained from an SEM image of a cluster of commer- 
cial 50 nm colloidal gold spheres. The gray scale values 
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were converted to a sample thickness varying between 
and 50 nm, and we assigned the complex index of refrac- 
tion of a 750 eV x ray photon going through an organic 
compound (PMMA). Here the numerical tests are done 
assuming periodic boundary conditions. These boundary 
conditions ensure that every region of the object ip is illu- 
minated with an equal number of overlapping frames, in 
other words the null space of Q is empty. We use frames 
width of 16x16, probe width=8, step size =5, number of 
frames=8x8... 64x64, RAAR algorithm, f3 = .75. 

The error metrics SF^e q used to monitor progress are: 



\a\\e F (*W) = \\[P F -I\ 



(19) 
(20) 



where I is the identity operator. This has to be compared 
to £q 7 the error w.r.t the known solution: 



\a\\e (zW) 



mm 



(21) 



where ip is an arbitrary global phase factor. 
We report the following observations 



# frames 


clock time(s) 


iteration 


el 


standard 








4x4 


0.7 


121 


<le-ll 


8x8 


1.4 


125 


<le-ll 


16x16 


4.9 


144 


<le-ll 


24x24 


26.3 


400 


4.3e-10 


32x32 


36.3 


400 


4.3e-4 


48x48 


90.7 


400 


3.4e-4 


64x64 


137.5 


400 


5.3e-3 


augmented 








4x4 


1.9 


138 


<le-ll 


8x8 


2.7 


141 


<le-ll 


16x16 


6.5 


138 


<le-ll 


24x24 


14 


134 


<le-ll 


32x32 


25.6 


139 


<le-ll 


48x48 


60.4 


142 


<le-ll 


64x64 


96.2 


149 


<le-ll 



TABLE II: Performance of projection algorithms using 
matlab R2012a 64-bit (maci64) (lapack version 3.3.1, 
MKL 10.3.5) on 2x2.2GHz Quad-core Intel xeon using 
frames of dimension 16 x 16 . 



• Fluctuating intensities: (Fig. 4), (Intensity 
fluctuation in this tests are 20%) By solving 
the new LSQ problem we obtain accelerated con- 
vergence and exact reconstructions every time we 
tested the problem, see (Fig. 4). No degradation 
(above numerical precision) introduced by intensi- 
ties perturbed by 20%. 

• Scaling: (Fig. 5) Improved convergence with 
larger problems. As we increase the number of 
frames, convergence slows down for standard pro- 
jection operators (parameters: m=16, Dx=4, k 
varies, n=k*Dx+m) 

• Multiplexing: (Fig. 5) Deconvolution of the in- 
coherent sum of frames translated by 3 times the 
illuminating beam width. 

• Calibration: (Fig. 6). Deconvolution of the in- 
coherent sum of frames translated by 3 times the il- 
luminating beam width, with unknown amplitude. 

• Positions: Fig. 7 Recovery of the positions per- 
turbed by an unknown factor randomly distributed 
between ±2.5 pixels. 

• Background: (Fig. 9), (||z (i) ||>*/||fe|| = 0.5. In 
Figure 8,9 we obtain exact reconstruction of the ob- 
ject and background (Background ratio ||a||/||b|| = 
1(T 6 . Exact recovery (within numerical precision) 
was obtained with reduced step size Sx = 3r. 
No degradation (above numerical precision) intro- 
duced by the background, nearly no influence on 
convergence rate. 



The dominant computational effort at experimentally 
practical size scale (k <C n) are the computation of frame- 
frame dot products and the linear solution of H\l. Each 
matrix entry Hij requires computation of the norm or 
dot product of a pair of overlapping frames, each mul- 
tiplied by the conjugate of the probe, normalized by 
(Q*Q) _1 and translated by — x^-). Since every frame 
has a fixed number of neighbors, this computation grows 
linearly with the number of frames k. Inversion of the 
(k x k) sparse matrix H, or the simple computation of 
H~ x l is typically a smaller computational problem than 
the computation of its entries. 



CONCLUSIONS 

The high redundancy in ptychographic data enables 
not only robust phase recovery but the recovery other 
parameters as well, such as the illuminating function 
itself 17 , position 13 ' 18 , coherence function 16 . Here we in- 
troduced a modified projection operator for the ptycho- 
graphic reconstruction problem that accounts for fluctu- 
ating intensities, position errors, or poor calibration using 
multiplexing illumination, and an unknown offset (back- 
ground) different for every pixel but constant throughout 
the acquisition (or vice versa). 

By optimizing the phase among each frame, we solve 
the phase problem at a resolution given by the step size 
between frames. This intra- frame phase optimization 
may be applied to merge subregions reconstructed in- 
dependently by distributed computer systems; For three 
dimensional objects, to merge two dimensional views re- 
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constructed independently into one three dimensional ob- 
ject. Finally, this intra- frame optimization could be ap- 
plied to multi-scale reconstructions where frames are di- 
vided in regions of Fourier space. 

The high redundancy in the data enables the identi- 
fication of background fluctuations over 10 6 x the signal 
and possibly more. We considered background fluctua- 
tions that are constant throughout a reduced dimension: 
either a flat offset that fluctuates for each snapshot, or 
an offset different for every pixel but constant through- 
out the acquisition. We introduced a parameter r] that 
scales intensities along the reduced dimension in the op- 
timization process that we found critical to achieve con- 
vergence. 

Our implementation of the background noise subprob- 
lem optimization has negligible effect on clock time or 
memory footprint. The subproblem associated with fluc- 
tuating intensities increases the clock time per iteration 
but improves the convergence rate for large scale prob- 
lems. 

More work is needed to include other sources of noise, 
and to establish the optimal frequency of communication 
and the amount of overlap between sub-reconstruction 
regions. 
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Appendix A: Solution to the augmented problem 



Here we set Yli=i c (i) = ^ anc ^ s °l ve 



min £(/0,c, A), where C = V \\Q^ - c (i) z (i) || 2 + A Vc (t) - k , 



(Al) 



where A is a Lagrange multiplier. We could also set c\ = 1 or any arbitrary number but numerical tests indicated 
improved performance of the former constraint. As we can see below, we can eliminate ip in by block factorization 
and solve: 



argmin - P Q )Diag (c) z|| 2 + A (^c, - fc) . 
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We note that G C are complex. To find the coefficients , we use the normal equation associated with the LSQ 
problem (Eq. Al) : 



( T,iQ*i)Q(i) -Qi*i -Q* k *k o\ 

-z\Q x z\z x ... 1 
: : : 



-zlQk 




We can partition the equation above as 



A B* 
B D 





ztz k 1 



(: 



where 



i=\ 



( -zlQi \ 



-zlQk 

V / 





1 


0/ 






( 


:) 


)■ 


(° 










( 













f4\ 



Ck 

V A/ 








D 



l\ 







... z* k z k 1 



V 1 



Using the block factorization 



I \ A 



B* 



BA~ l I/O D — BA~ l B* 



We obtain: 



(D — BA~ l B*) 



c 



/ v> \ 



1 0/ 



/ o 



We can express the Schur complement: 



D — BA~ l B* 



H 1 

1* 



/1\ 



\1/ 



\c k J 



(A2) 



where 



HiJ — Z (i) 



(i) 



V 



(0 



/ 



Note: Hc = V c ||(/ - P Q )z Diag (c) f 

By block-wise inversion, of Eq. A2 we obtain c and the scaling factor A from a sparse linear equation: 

Hc=-Xl, -A= (l*^- 1 !)- 1 ^ 
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Appendix B: Taylor expansion 



We consider the case in which the probe w is translated from the input coordinate by an unknown distance £. 
We call the unknown illumination matrix used to generate the data. To determine the illumination matrix, we 
determine the parameter £ so that the error sq^ is minimized: 

i i|2 



arg min 1 1 [/ • 
£ 11 



(Bl) 



We restrict initially to 1-dimensional translations. Given the illumination function w, we can compute the first and 
second order derivatives with respect to translation. We denote by Qi,Ri, Si the illumination matrices that extract a 
frame out of an image and multiplies by (w(x),d x w(x), d%w(x)) respectively. Q,R, S are tall and skinny matrices of 
the same size as Q discussed earlier, with identical location of the non-zero entries. The probe is perturbed to second 
order as follows: 

Q^ = Q + £R + £ 2 S. 

where £ is diagonal and real, and represents the translation distance different for every frame i. We minimize 

-l 



arg mm 



By Taylor expansion: 



/-(Q + £R + £ 2 S) (Q + £R + £ 2 S) (Q + £R + £ 2 S) (Q + £R + £2 S ) 



(B2) 



1 



Q*Q 



(R*£Q + Q*£R + 0(£ 2 )) 



l 



Q*Q 



the term 0(£ 2 ) includes other second order terms that we will not need. We write the expansion of the residual in 
Eq. (Bl)/o + /i(0 + / 2 (£ 2 )as: 

fo = [I- Pq}z = z (B3) 
We define <p* = qtqQ*, and express the first order as 

h = \~W - 0R*£ + <MR*£Q + Q*£R)<t>*\ z 
by defining z R = R^*z, using Pq = Q(p* = and rearranging: 

fi = -?z R - </>R*£(z - Q<f z) + 0Q*£z R 

= -(1 - P Q )£z R - 0R*£z (B4) 
using the equality z(J — Pq) = z, setting O r = R</>*</>R* = RqIqR* and rearranging we obtain: 



h = - [CRq^qR*C + C 2 S0 + QO (£ 2 ) 
= -£0 R £z - C 2 zs + QO (C 2 ) z 
where zs = S0*z. We rewrite Eq. (Bl) above as: 

|2 



(B5) 



II [/ - p Qi \ zf = fsfo + fsfi + /r/o + /r/i + /o/ 2 + /i/o + o(e) 

we note that z*Q = z*0 = 0, set zq = Q^z and obtain the first and second order terms of Eq. (Bl): 

/o/i + /r/o = -z*£z R , - Zr £z (B6) 
/r/i + /o h + / 2 7o = z R C(/ - Pq)?z r + z£O r £z - z*£O r £z _ ? * ^ Zs _ z * ?0r ^ _ 

= z R £ 2 z R - z*£ 2 z s - z^ 2 z - z*£0r£z + z Q £0 R £z Q - z r £Pq£z r , 
By using the definition of zq, z r , Pq and R it is easy to show that Zq£O r £zq = z r £Pq£z r and simplify as: 

fih + /0V2 + / 2 */o = z^ 2 z R - z*£ 2 z s - z^ 2 z - z*^z, (B7) 
~ z^ 2 z R - z^Or^z. (B8) 

We note that — z*£ 2 zs — Zg£ 2 z is a second order correction if we assume that z is in that range of an unknown 
for small £. By setting d^.\\ • || 2 = 0, we obtain the linear equation: 

( 2 ( Z R* Z R* ~ zJ z s, ~ z k*i) S ij ~ Z i°^ij z j ~ tfOnjiZi) = zJz R . + z^ 
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FIG. 2: object ip used to simulate diffraction data) 



FIG. 3: Absolute value of the probe \w(r)\ used in 
simulations (16x16 pixels) 
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FIG. 4: Convergence rate with an 10-error of ±20%. (left) old projection operator (right) new projection operator, 
(bottom) reconstruction from data with 10-error, and solution (reconstruction using the new projection operator is 

within the computer numerical precision). 
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errors in ptychography. Ultramicroscopy, 120(0) :64 - 72, 2012. 
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(a)8x 8 frames (b) 16x16 frames (c)32x 32 frames (d)64x 64 frames 
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(e)8x8 frames (f)16xl6 frames (g)32x32 frames (h)64x64 frames 



FIG. 5: Convergence rate (ep, £q, £o vs number of iteration i) for (top) regular reconstruction, (bottom) using 

augmented projection (m = 16 and step size x\ — X2 = 3) 
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FIG. 6: Convergence rate with incoherent illumination of 4 beams, separated by 3x the probe width (FWHM) using 
standard projection algorithms (top left), with intermediate phase optimization (top-right), phase and 
amplitude (bottom-left), and phase and amplitude with initial amplitude error of 20% (bottom right), frame width 
16 x 16, 16 x 16 frames, step size 3.5 pixels close packing with ±1 pixel known random perturbations. 
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FIG. 7: Reconstruction with position errors using the method described in section V, where e^ = ||£ — £o||/||£o||> 

and the perturbations in position are randomly distributed with (£o) = \ J2i II& II = 2.5 resolution elements, 
(number of frames: 16x16, frame dimensions 32 x 32, step size: 3.5 pixels, hexagonal packing with known random 
perturbations of ±1 pixels and unkown £ random perturbations ). 




FIG. 8: Two measured intensities with additive background (SNR=0.5). In a separate test the diffraction data was 
buried by the background ( (in other figures not included background was 10 6 x the signal). 
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FIG. 9: (top) reconstructed image with background optimization (left) and without (right), (bottom) reconstructed 

background (left), convergence behavior (right) 



